Tunneling of anyonic Majorana excitations in topological superconductors 
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We consider topological superconductors and topological insulator/superconductor structures in 
the presence of multiple static vortices that host Majorana modes and focus on the Majorana tun- 
neling processes between vortices. It is shown that these tunnelings generally lift the degeneracy of 
the many-body ground state in a non-universal way, sensitive to microscopic details at the smallest 
length-scales determined by the underlying physical problem. We also discuss an explicit realiza- 
tion of the Jackiw- Rossi zero- mode in a topological insulator/superconductor structure with zero 
chemical potential. In this case, the exact degeneracy of the many-anyon ground state is protected 
by an additional chiral symmetry and can be linked to the rigorous index theorem. However, the 
existence of a non-zero chemical potential, as expected in realistic solid state structures, breaks 
chiral symmetry and removes protection, which leads to the degeneracy being lifted. Finally, we 
discuss the implications of our results for the collective states of many-anyon systems. We argue 
that quantum dynamics of vortices in realistic systems is generally important and may give rise 
to effective time-dependent gauge factors that enter interaction terms between Majorana modes in 
many-anyon systems. 
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I. INTRODUCTION 

Topological quantum computation hinges on the ex- 
istence of non-Abelian excitations, which arise in cer- 
tain topological phases of matter pj. The first exam- 
ple of such a state is Fractional Quantum Hall (FQH) 
state with v = 5/2 filling fraction. This state is be- 
lieved to be described by the Moore-Read Pfaffian wave 
function Q which supports non-Abelian excitations §■ 
It was shown later that the Moore-Read Pfaffian wave 
function of composite fermions is related to the BCS wave 
function with p-wave pairing establishing the con- 

nection between FQH state v — 5/2 and p x + ip y super- 
conductors. Certain vortex excitations in chiral p-wave 
superconductors carry zero energy modes and obey non- 
Abelian statistics [1, 0, HI • These topologically protected 
zero modes can be occupied by Majorana fermions and 
are responsible for topological ground state degeneracy. 
Namely, 2n vortices with zero modes residing in vortex 
cores span 2"~ ^dimensional Hilbert space. Non-Abelian 
statistics of vortices can also be derived within this frame- 
work [8| . Making use of their intrinsic non-local quantum 
entanglement, vortices carrying Majorana modes can be 
exploited to realize topological qubits which are inher- 
ently decoherence-free and are protected against smooth 
local perturbations, thus providing a very appealing plat- 
form for fault-tolerant topological quantum computa- 
tion 0-tni- 

From the perspective of experimental realization of 
topological phases, there is some preliminary evidence 
that v = 5/2 FQH state may have non-Abelian exci- 
tations [12l4l4| . Spin-triplet p x + ip y pairing superfluid- 
ity/superconductivity are believed to occur in A-ph&se of 
superfluid 3 He [HI, [l6[ and strontium ruthenates fl7l - [2~fl ] 



in which half-quantum vortices are non-Abelian [8|, l22j . 
There are also proposals to realize chiral p-wave super- 
fluids in ultracold atom systems [23U281 ] . The theoreti- 
cal description of these systems essentially falls into the 
category of a spinlcss p x + ip v superconductor. Apart 
from these examples, there also exist a number of pro- 
posals involving various heterostructures of a three di- 
mensional to polo gical insulator (TI) and a supercon- 
ductor (SC) [29l |30|, semiconductor and superconduc- 
tor [3l|, [32J and superconductor and ferromagnet [H, [34j . 
These systems seem to be more experimental accessible. 
We also note that there are proposals to realize Majo- 
rana fermions in one-dimensional systems (35l - [39| and on 
the surface of a three dimensional Z2 topological super- 
fluid/superconductor [4CH43 | . 

We note here that the emergent Majorana excitations 
in these physical systems are closely related to zero- 
modes that have been long known in the context of 
high-energy physics, where chiral fermions in the pres- 
ence of topological defects (domain walls, vortices, etc.) 
give rise to massless excitations within the topological 
defects [44J. One such example is a Jackiw- Rossi zero- 
mode that was predicted to appear in a vortex-chiral- 
fcrmion system. An important ingredient of the original 
Jackiw-Rossi model 45] is the existence of a conserved 
chiral current, which allows one to enumerate zero-energy 
modes according to their chirality. It was subsequently 
established by Weinberg [46| that for the Dirac opera- 
tor the difference between the number of zero modes of 
opposite chirality is given by the winding number of the 
superconducting order parameter phase. Thus, the con- 
servation of the chiral current enables the use of a pow- 
erful Atiyah-Singer-type index theorem [U EH, which 
relates a number of zero modes to the total topological 
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charge of the vortex configuration. In particular, the the- 
orem ensures that the exact degeneracy of the many-body 
ground state is preserved and no tunneling process can 
possibly lift it. As discussed in Refs. [29|, [Jjj, a Hamil- 
tonian considered by Jackiw-Rossi can be realized in a 
Tl/supcrconductor heterojunction. However, as shown 
below a realistic solid-state structure of this type is gen- 
erally described by a low-energy theory, which has the 
exact chiral symmetry only if the chemical potential of 
excitations is exactly zero. Due to the non-trivial way 
the chemical potential enters the BdG equations, there 
is no conserved chiral current, and one can not enumer- 
ate zero modes by the their chirality anymore. Therefore, 
the connection between the analytical and topological in- 
dex established in Ref. [!(| does not apply in this case, 
and intervortex tunneling leads to lifting of the ground- 
state degeneracy of a many-anyon system for any finite 
chemical potential. Clarification of the applicability of 
the index theorem to the non-relativistic topological su- 
perconductors is the main result of the paper. 

Understanding the fate of ground-state degeneracy of 
many-anyon system in realistic solid-state structures is 
a difficult problem of fundamental importance and of 
relevance to practical realization of topological quantum 
computing. In this paper we address one mechanism that 
may lift the ground state degeneracy associated with the 
tunneling processes between spatially separated vortices. 
The presence of the bulk gap protects ground state de- 
generacy from thermal fluctuations at low temperature 
leaving out only processes of Majorana fermion quantum 
tunneling between vortices. Generic features of tunneling 
of topological charges have been explored recently [50j . 
The lifting of ground state degeneracy due to intervor- 
tex tunneling for a pair of vortices have been studied 
numerically for v = 5/2 quantum Hall state [HI, [52j, 
Px + iPy superconductor [53| and Kitaev's honeycomb lat- 
tice model jUj]. Analytical calculation has been carried 
out for the model of spinless p x +ip y superconductors [55$ . 

Generally energy splitting due to intervortex tunneling 
is determined by the wave function overlap of localized 
Majorana bound states. In this paper we calculate the 
splitting for both spinless p x + ip y superconductor and 
a model of Dirac fermions interacting with the scalar 
superconducting pairing potential realized in a TI/SC 
heterostructure. In both cases, besides the expected ex- 
ponential decay behavior, it is found that the prefactor 
exhibits an oscillatory behavior with the intervortex dis- 
tance which originates from the interference of different 
bound state wave functions oscillating with the Fermi 
wave length. This is generic situation for weak coupling 
superconductors where the Fermi energy Ep is much 
larger than the superconducting gap A. In this paper, we 
also consider several cases where the Fermi wavelength is 
much larger than the coherence length. This scenario is 
relevant, for example, for TI/SC heterostructure as well 
as some other systems involving the proximity-induced 
superconductivity. When chemical potential is tuned to 
the Dirac point (/i — > 0), we find indeed that the split- 



ting in TI/SC heterostructure vanishes. This fact can 
be attributed to an additional symmetry possessed by 
the system at = - the chiral symmetry as discussed 
above. 

The paper is organized as follows. In Sec. |H] we re- 
view the Bogoliubov-de Gennes equations for spinless 
Px + ipy superconductor as well as TI/SC heterostruc- 
ture and show that there are Majorana zero energy so- 
lutions localized at the vortex core. Then using these 
bound state wave functions, we calculate energy split- 
ting of zero-energy states due to intervortex tunneling in 
Sec. IIIII We present our main results in Sec. II VI by inter- 
preting the explicit splitting calculations presented in the 
previous section from the perspective of the index theo- 
rem which establishes the relation between zero-energy 
modes and topological index of the order parameter. Im- 
plications of our result for topological quantum compu- 
tation and interacting many anyons system are discussed 
in Sec. El 



II. ZERO ENERGY BOUND STATES IN 
SUPERCONDUCTING VORTEX CORES 

In this Section, we review the analytic solution of the 
Bogoliubov-de Gennes equation for a zero-energy Majo- 
rana mode in a p + ip superconductor and at a topologi- 
cal insulator/superconductor interface. The main results 
here are Eqs. © and (|2"5j) and they are used further in 
Sec. Ill to calculate the energy-level splitting due to tun- 
neling. 



A. Bound states in p x + ip y superconductors 

We start with the mean-field Hamiltonian for spinless 
Px + ipy superconductor 

Hbcs = J d 2 r^( r ) _ ^ ^( r ) + 



d 2 rdV 



(1) 



^ t (r)A(r,r')V' t (r / ) +h.c 

where the gap operator A(r,r') is given by (5|| 

A(r, r') = -lA (l±£^ (d x , + id y ,)S(r - r'). (2) 

This Hamiltonian can be derived as a continuum limit 
of a lattice model of spinless fermions 57]. To diag- 
onalize this Hamiltonian we perform Bogoliubov trans- 
form -0(r) = J2n [% u n(r) + jX v n( T )] where n labels 
different quasiparticle eigenstates. Canonical commuta- 
tion relation ["HBCS,7n] = —E n j n yields corresponding 
Bogoliubov-de Gennes(BdG) equation 



'HBdG 



u(r) 
v(r) 



E 



u(r) 
v(r) J ' 



(3) 
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where BdG Hamiltonian reads 



HBdG- 



2m 



■fj, 



A(r), d x + id v } 



1 



\— pA^r),^-^} 



2m 



(4) 



with anti-commutator being defined as {a, b} = (ab + 
ba)/2. 

Before discussing the zero-energy solutions of BdG 
equations, it is instructive to review the symme- 
tries of the Bogoliubov-de-Gennes Hamiltonian (j4}. 
Particle-hole symmetry of BdG Hamiltonian follows from 
{Sj^Bdc} — where 3 = t x K with K being com- 
plex conjugation operator and t x being Pauli matrix in 
Nambu(particle-hole) space [58]. Besides particle-hole 
symmetry BdG Hamiltonian (Q} has no other generic 
symmetries. Thus, it is a typical example of the symme- 
try class D in the general classification scheme of topolog- 
ical insulators and superconductors [H, |6(| • The direct 
consequence of particle-hole symmetry is that the eigen- 
states of -ffBdcQ come in pairs, with opposite eigenen- 
ergies. That is, if ^ = {u n ,v n ) T is a solution of Eq. (O 
with eigenvalue E n , then E^f = (w*,m*) t must be a so- 
lution with the eigenvalue (—E n ). Particularly, a non- 
degenerate zero energy state must obey the following 
constraint: 3^ = X^. Because S 2 = 1 which im- 
plies |A| = 1, A must be a pure phase A = e l6 . Wc 
can make a global gauge transformation and introduce 
Nambu spinors as = exjp(-i9/2)^i then 3*' = 
Thus, a non-degenerate zero energy state should always 
satisfy u* — v. The corresponding quasiparticle operator 



7' 



d 2 r [u(r)^ t (r) + v(r)tp(r)] 



is then self-Hermitian obeying 7 = -yt, i.e. 7 is a Majo- 
rana fermion operator. 



We will now show that such zero energy states appear 
in the cores of vortices in chiral p-wave superconduc- 
tors. The localized states in the vortex cores are known 
as Caroli-de-Gennes-Matricon states (CdGM) [6l|. In 
conventional s-wave superconductors all CdGM states 
have non-zero energies Q. However, due to the chi- 
rality of the order parameter p x 4- ijh, superconductors 
can host zero-energy bound states p], HH, [HI, [H, l63j |. 
The non-degenerate zero-energy bound states are topo- 
logically protected by the particle-hole symmetry. The 
existence of the zero energy solution in the vortices of 
the chiral p-wave superconductors can be demonstrated 
explicitly by solving BdG e quat ions [HI]. Similar to the 
s-wave superconductors [Hll |64|. a vortex with vorticity 
l(i.e. I flux quanta hc/2e is trapped) can be modeled as 
A(r) = f(r)e llip where (p is the phase of the order pa- 
rameter and f(r) is the vortex profile which can be well 
approximated by f(r) = A tanh(r/£) [13] ■ Here A is 
the mean-field value of superconducting order parameter 
and £ = vp/Ao is coherence length. We will mainly focus 
on the case 1 = 1. Taking advantage of rotational sym- 
metry, BdG equation can be decoupled into angular mo- 
mentum channels. The wave function can be written as 
# TO (r) = e imv (e lv u m (r),e- iip v m (r)). As argued above, 
a non-degenerate zero mode requires 3^ cx 'J which can 
only be satisfied for m = 0. The radial part of the BdG 
equations in m — channel then reads 



"2^ + ^-4,) 

■£/(r) \(dr + ±) + 



f'(r) 
2 



2m 




= 0. 



(5) 



Given that the radial part of the BdG equation is 
real, one can choose uo(r) and vo(r) to be real. Then 
the condition 3^o = ^0 reduces to vq — Auo with A ~ 
±1. Using this constraint, the differential equation for 
Mo becomes: 



1* • 



— 2mjx- 



2A 

vf 



fid, 



1 

2^ 



L 

2 



M =0. 



One can seek the solution of the above equation in the 
form u(r) = x(r) exp [A J Q dr' f(r')\ leading to 



X "+- +(2m^-C-\ 



x = o. 



(6) 



Here the profile f(r) — Aotanh(r/£) vanishes at the ori- 



gin and reaches Ao away from vortex core region. For 
our purpose, it's sufficient to consider the behavior of so- 
lution outside the core region where f(r) is equal to its 
asymptotic bulk value Ao . It is obvious now that A = — 1 
yields the only normalizable solution. 

When Ao < 2m[iv F which is the case for weak- 
coupling BCS superconductors, Eq.® becomes first or- 
der Bessel equation. Thus, the solution is given by Bessel 
function of the first kind J n (x): 



X(r) = M 1 J 1 {r^2m t x-A 2 /v 2 F ), 



(7) 



where Afi is the normalization constant determined by 
the following equation 4ir J rdr |Mo(r)| 2 = 1. In the op- 
posite limit Aq > 2mfj,Vp, the solution of Eq. © is given 
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by first order imaginary Bessel function: 



x (r) =Af 2 I 1 (r^A 2 /v 2 F -2mii) 



(8) 



The function I n (r) diverges when r — >■ oo. But the 
radial wave function Uo(r) remains bounded as long as 
fi > 0. This is consistent with the fact that \x = sepa- 
rates Abelian topological phase (/i < 0) and non-Abelian 
phase(/i > 0) |5[. The critical value Aq = 2mnv F cor- 
responds to closing of the bulk gap. At this point the 
notion of a localized bound state becomes meaningless. 
Indeed, at this point the solution of Eq. © scales as 
X(r) oc r. 

We summarize this section by providing an explicit 
expression for zero energy eigenfunction: 



*o(r) = x(r)exp 



vfJ 



dr'/(r') 



(9) 



where x( r ) is given by Eq.([7]) for Aq < 2mfiv F and Eq.® 
for Aq > 2m[iv F . 

Using the zero energy solution obtained for one vortex 
one can be easily write down wave function for multi- 
ple vortices spatially separated so that tunneling effects 
can be ignored. Assume there are 2N vortices pinned at 
positions Rj ,i = 1, . . . , 27V. The superconducting order 
parameter can be represented as 



:2.\ 



A(r) = n/(r-Ri)exp[i5>i(r) 



(10) 



where (pi(r) = arg(r — Rj). Near the k-th vortex core, 
the phase of the order parameter is well approximated 
by y>fc(r) + Slfc with Slj, = J2i^k <fii(R-k) which is accurate 
in the limit of large intervortex seperation. Then in the 
vicinity of k-th vortex core, a zero energy bound state 
can be found 15611: 



r x(n=)exp 



vf Jo 



dr'/(r') 



x exp 



£1^ 

i\ fk + — I T > 



(11) 



where — \r — Rfe|. Correspondingly, there are 2N Ma- 
jorana fermion modes localized in the vortex cores. They 
can be combined pairwise to form N Dirac fermions. 
Specifically, two Majorana fermions 7 and 7j localized 
in i-th and j-th vortex cores, respectively, are combined 
into a Dirac fermion: 



1 



V2 



(12) 



These N Dirac fermions can be occupied or unoccu- 
pied, allowing for enumeration of all degenerate ground 
states [H. 



B. Bound states in the Dirac fermion model 
coupled with s-wave superconducting scalar held. 

We now discuss the zero energy bound states emerging 
in the model of Dirac fermions interacting with the super- 
conducting pairing potential. This model is realized at 
the interface of a 3D strong topological insulator having 
an odd number of Dirac cones per surface and an s-wave 
superconductor [29j. Due to the proximity effect an in- 
teresting topological state is formed at the 2D interface 
between the insulator and superconductor. We will now 
discuss the emergence of Majorana zero energy states at 
the TI/SC heterostructure [29]]. This model was also con- 
sidered earlier in the high-energy context by Jackiw and 
Rossi f45] . 

Three dimensional time-reversal invariant topological 
insulators are characterized by an odd number of Dirac 
cones enclosed by Fermi surface [65T - lo7| . The metallic 
surface state is described by the Dirac Hamiltonian. The 
non-trivial Z2 topological invariant ensures the stabil- 
ity of metallic surface states against perturbations which 
preserve time-reversal symmetry. When chemical poten- 
tial fx is close to the Dir ac p oint the TI/SC heterostruc- 
ture can be modeled as [29l lo8j: 

n = ^(va-p- + A^j+h.c, (I 3 ) 

where ijj = (ip-f, ipi) T and v is the Fermi velocity at Dirac 
point. The Bogoliubov-de Gennes equations are given 
by: 



%BdG 



H BdG *(r) = £tf (r) 
a ■ p — /i A 

A* -<7-p + fl 



(14) 
(15) 



where ^(r) is the Nambu spinor defined as ^ = 
(u^, uj,, — Vf) T . At /i = the BdG Hamiltonian above 
can be conveniently written in terms of the Dirac matri- 
ces: 



H B dG = (laPa + T a n a ) . 



(16) 



a=l,2 



Here 7 a and T a are 4x4 Dirac matrices defined as 
7i = 72 = <T v T Zl and Ti = T X ,T 2 = r y and n = 

(ReA, — ImA). One can check that these matrices sat- 
isfy the following properties: {7 a ,7fc} = {r^T,,} = 5 ab 
and {r a ,7(,} = 0. The fifth Dirac matrix 75 is given by 
75 = -7i72rir 2 = T z a z . 

As in the case of spinless p x + ip y case, we first discuss 
the symmetries of Eq. (|15[) . The particle- hole symmetry 
is now 3 = UyTyK where r are Pauli matrices operating 
in Nambu (particle-hole) space. The difference with the 
previous case is the presence of time-reversal symmetry: 
6 = i<jyK, [9,"HBdG] = in this model. Moreover, when 
fi = there is additional chiral symmetry in the model 
which can be expressed as {7 5 ,"HBdG} = 0. We will 
see that the chiral symmetry has important implications 



5 



for degeneracy splitting. Bogoliubov quasiparticles are 
defined from solutions of BdG equations as 



r)W(r). (17) 



If we require 7 to be a Majorana fermion, i.e. 7 = 7^, 
the necessary and sufficient condition is v a = u* up to a 
global phase. 

A vortex with vorticity / can be introduced in the or- 
der parameter as A(r) = f(r)e llip . Rotational symmetry 
allows decomposition of solutions into different angular 
momentum channels: 



*m(r) 



e-**/4 Xt ( r ) \ 



e- i7r / 4 774.(r)e-^ 



i7r/4 



?7 t (r) 



We define $0 = (xt) X4.) *74.i ^t) 7 f° r convenience. 

Similar to the previous analysis, we first look for non- 
degenerate Majorana zero-energy state. The Majorana 
condition 5^ oc ^ fixes the value of to to be for odd 
I. For even I, there is no integer to satisfying Majorana 
condition so no Majorana zero mode exists. The radial 
part of BdG equation then becomes 



n r f(r) 
f(r) 



Ur = 



(3 yl~Lf O y 



*o(r) = 

V (d r + 



(18) 
(19) 



Here "I/o is assumed real. Since we are interested in non- 
degenerate solution, must be simultaneously an eigen- 
states of (T y Ty (particle- hole symmetry). This condition 
implies that 771 = — Xxt^Vl — Xxi where A = ±1. Tak- 
ing into account above constraints 4x4 BdG equation 
reduces to 



-H V (d r + ^)+ Xf 

-V (d r -f)~ Xf -/I 



= 0. 
(20) 

The solution of the above equation can be easily obtained 
for 

ty-Htt&y''"""- (21) 

Obviously, we should take A = 1 to make radial wave 
functions normalizable. Here A3 is the normalization 
constant whose analytical form is given in (|A7[) . 

The case of \x = is special due to the presence of an 
additional symmetry of BdG Hamiltonian - chiral sym- 
metry. Imposing the boundary condition at r — > that 
wave function must remain finite, for I > the solution 
of Eg. ([20]) becomes 



So dr ' f(r') 



(22) 



and if I < 0, m + 1 < 0, 







.-(m+l) 



■Jo r <lr'/(r') 



(23) 



Particle-hole symmetry combined with the analyticity of 
the above solutions at r constraints the integer 
< to < I — 1 . Similar conditions apply to the negative I 
solutions. Thus, there are exactly \l\ zero-energy modes 
for any I as found previously in Ref. [45J] . If Hs even, all 
these solutions are Dirac fermionic modes. However, if I 
is odd, there are I — 1 Dirac fermionic modes and one Ma- 
jorana zero-energy mode as given in Eqs. (1221) and (|23p . 
Because the chiral symmetry also relates eigenstates with 
positive energies to those with negative energies which 
follows from 7 5- Hbc1G7 5 = — 'HBdGj one can always re- 
quire the zero-energy eigenstates to be eigenstates of 7 5 . 
The wave function in Eq. (|2"2"| is an eigenstate of 7 s with 
eigenvalue I while wave function (f2"3"f has eigenvalue —1. 
We define eigenstates of 7 5 with eigenvalue ±1 as ± chi- 
rality. 

To summarize we have obtained the Majorana zero- 
energy bound state attached to the vortex with odd vor- 
ticity: 



*o(r) 



_ e i(J-l)v»/2 



/ 



V- 



t/4 



e W4 
i7r/4 



Xf( r ) \ 



Jir/4 



Xi(r)t 
Xl{r)e 



X\(r)e 



—Hip 
-i(l-l)tp 



(24) 



Generalization to the case of many vortices is straight- 
forward. Order parameter with 2N vortices pinned at 
Ri is already given in (fTU| . Assuming that they are well 
separated from each other, we can find an approximate 
zero-energy bound state localized in each vortex core: 



e-™' A x^n) \ 

e" /4 Xi(n)e^ 



-iir/4 



-ilipi 



y- e W4 Xt ( r .) e -i0-i) Vi y 



(25) 

the construction of N Dirac fermions and 2 N 1 ground 
state Hilbcrt space are the same as the case of spinless 
Px + iPy superconductors. 



III. DEGENERACY SPLITTING DUE TO 
INTERVORTEX TUNNELING 

The ground state degeneracy, which is crucial for topo- 
logical quantum computation with non-Abelian anyons, 
heavily relies on the assumption that intervortex tunnel- 
ing is negligible. When tunneling effects are taken into 
account zero energy bound states are usually splitted and 
the ground state degeneracy is lifted. Besides, the sign 
of energy splitting is important for understanding many- 
body collective states [69| . 

We now discuss a general formalism to calculate the 
energy splitting. We focus on the case of two classical 
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vortices each with vorticity I = 1 located at certain fixed 
positions Hi and R 2 . To develop a physical intuition, it 
is useful to view a vortex as a potential well, which may 
host bound states including zero-energy states, while the 
regions where superconducting gap is finite play the role 
of a potential barrier. Therefore, the two- vortex problem 
resembles the double-well potential problem in single- 
particle quantum mechanics (sometimes referred to as 
the Lifshitz problem in the literature [z3|). The solution 
to this simple problem in one-dimensional quantum me- 
chanics is readily obtained [70| by considering symmet- 
ric and antisymmetric combinations of single- well wave- 
functions (which can be taken within the quasiclassical 
approximation for high barriers) and the overlap of these 
wave-functions always selects the symmetric state as the 
ground state in accordance with the elementary oscilla- 
tion theorem (i.e., the ground state has no nodes). We 
note that both quasiclassical approximation and the Lif- 
shitz method are not specific to the Schrodinger equation, 
but actually represent general mathematical methods of 
solving differential equations of certain types. Moreover, 
these methods can be applied to rather generic matrix 
differential operators, and such a generalization has been 
carried out by one of the authors in a c omp letely differ- 
ent context of magnetohydrodynamics, |7l| where inter- 
estingly the relevant differential operator appears to be 
mathematically similar to the BdG Hamiltonian. These 
considerations suggest that one can use the generalized 
Lifshitz method to obtain the splitting of zero modes of 
the BdG equations, by considering certain linear com- 
binations of the individual Majorana modes in the two 
vortices and calculating their overlap, which reduces to 
a boundary integral along a path between the two vor- 
tices. Also, if the inter-vortex separation is large, one 
can use the semiclassical form of the Majorana wave- 
functions (effectively their large-distance asymptotes) to 
obtain quantitatively accurate results. Let us note here 
that apart from a technically more complicated calcula- 
tion that needs to be carried out for the BdG equation, 
another important difference between this problem and 
the simple Lifshitz problem is that we can not rely on 
any oscillation theorem and there is no way to determine 
a priori which state has a lower energy. As discussed be- 
low, this "uncertainty" is fundamental to this problem 
and is eventually responsible for a fast-oscillating energy 
splitting with intervortex separation. 

With the two zero-energy eigenstates \&i and ^2 lo- 
calized at Ri and R2 (given by Eq. (fTTj) for spinless 
Px + ipy superconductor and by Eq. (f!H|) for TI/SC het- 
erostructure), we can construct approximate eigenstate 
wave functions in the case of two vortices: "J ± = -^(\l/i± 

e la ^2) analogous to the symmetric and anti-symmetric 
wave functions in a double-well problem with energies 
E±, respectively. The phase factor e la can be deter- 
mined from particle-hole symmetry which requires that 
new eigenstates ^4. with energy E + = SE and 4"_ with 
energy E- = —SE be related by S\l/+ = Since 
^1 and ^2 are real (Majorana) eigenstates, one finds 
S* + = ^(^l + e~ ia ^ 2 ) = Thus, one arrives at 



e 2ia _ 2 which fixes a = ±tt/2. In the rest of the text 
we take a — n/2 for convenience. The corresponding 
quasiparticle operator can be identified with the Dirac 
fermion operator. We explicitly show this for the case of 
spinless p x + ip y superconductor: 



-172 ) 



d 2 r 



0- 



V2 



V2 



(26) 

Therefore c (&) annihilates (creates) a quasiparticle on 
energy level E+. The original two fold degeneracy be- 
tween state with no occupation c|0) = and occupied 



1 



|0) is lifted by energy splitting £4 



To calculate the energy of we employ the stan- 
dard method based on the wave function overlap [70| . 
Suppose the two vortices are placed symmetrically with 
respect to y axis: R x = (-R/2,0) and R 2 = (-.R/2,0). 
BdG equations are H B dG^+ = E + ^ + ,H B dG^i = 0. We 
then multiply the first equation by 4" J and second by 4^, 
substract corresponding terms, and integrate over region 
S which is the half plane x € (0, 00), y £ (—00, 00) arriv- 
ing finally at the following expression for E + : 



E, 



Jg d 2 r ^\n B dG^+ - Jg d 2 r ^^BdG^i 
/ s d2r*i* + 



(27) 



This is the general expression for the energy splitting 
which is used to evaluate E + in p x + ip y SC and TI/SC 
heterostructure. 



A. Splitting in spinless p x + ip y superconductor 

We now calculate splitting for two vortices in spinless 
Px + iPy superconductor. The denominator in Eq. (|27[) 
can be evaluated quite straightforwardly J s d 2 r 4 r i4 r + ~ 
With the help of Green's theorem the integral 
over half plane in the numerator can be transformed into 
a line integral along the boundary of S, namely the y 
axis at x = which we denote by 9S: 

E+ = mJ ^9 ( s ) 9' ( s ) cos 2tp 2 cos ^ 2 

.9 ( s ) ■ o • 9 is) 

H sin 2ip 2 sin ip 2 7 — 

s £ 



where s = y/ (R/2) 2 + y 2 , tany>2 = 2y/R. The function 
g(s) is defined as g(s) = x(s) exp(— s/£). 

First we consider the regime where Aq < 2m^iVp and 
radial wave function of Majorana bound state has the 
form ([7])- We are mainly interested in the behavior of 
energy splitting at large R £ with £ being the coher- 
ence length, where our tunneling approximation is valid. 
Another length scales in our problem is the length cor- 
responding to the bound state wave function oscillations 
k = y/2m[L — Ap/up. In the limit R ^> max(£; -1 ,£) 
upon evaluating the integral (|2"5)) we obtain 
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87V? / A 2 
7T m I 1 + A 2 



1/4 



: exp 



cos(fci? + a) — — sin(fci? 
A 



2(l + A 2 ) 1 /4 



(29) 



r 



where A = fc£, 2a = arctanA and A/"i is the normalization 
constant defined in Eq.([7]). The expression of 7V*i is given 
in Appendix |A1 and has the following asymptotes for A 3> 
1 and A < 1: 



/V 2 



A > 1 




(30) 



The exponential decay is expected due to the fact that 
Majorana bound states are localized in vortex core. In 
addition, the splitting energy E + oscillates with intervor- 
tex seperation R which can be traced back to interference 
between the wave functions of the two Majorana bound 
states since they both oscillates in space. 

Of particular importance is the sign of splitting as 
noted in Ref. [69]. It determines which state is ener- 
getically favored when tunneling interaction is present. 
If E + > 0, |0) is favored whereas E + < favors |1). We 
note here that the definition of states |0) and |1) relies 
on how we define the Dirac fermion operator c and c'. 
Due to the presence of a constant term together with 
trigonometric function, the sign of splitting can change. 
To figure out when the sign oscillates, we require the 
amplitude of the trigonometric part is greater than the 
constant part which gives 



4 2(1 + A 2 ) 1 / 4 
— > — 

A 2 A 



Solving this inequality yields A = k£ > 8. Therefore in 
this parameter regime the sign of splitting changes with 
distance R. Otherwise the splitting still shows oscillatory 
behavior but the sign is fixed to be positive. 

In weak-coupling superconductors where Ao < Ej? or 
equivalently kp£, 3> 1, the expression for the energy split- 
ting l|29p can be considerably simplified. In this case, 
(i ~ £p and k ~ hp- Keeping only terms that are leading 
order in (^fC) -1 , we find 



E, 



2 A cos(k F R+ f) / R 



(31) 



which is the expression reported in Ref. [551 ] . A 
similar expression for splitting of a pair of Majorana 
bound states on superconductor/2D topological insula- 
tor/magnet interface is found in Ref. [721 ] . 

Next we consider a different limit A 2 , > 2mp,Vp in 
which the wave function of Majorana bound state for a 
single vortex doesn't show any spatial oscillations. Thus, 
we expect that tunneling splitting will show just an expo- 
nential decay without any oscillations. The wave function 



(|8]) grows exponentially when r — > oo: 

1 



A-iir 



X (r) ~ -^e fc ° 
\ r 



with fco = v ^o/ v f — 2m/i. The overall radial wave 
function decays exponentially ~ exp(— k'r) where k' = 
l/£ — fco- In this case, the tunneling approximation is 
only valid for k'R ^> 1 since bound state wave function 
is localized approximately within distance 1/k' to vor- 
tex core. The resulting energy splitting monotonically 
decays: 



E, 



2A/f 



3 

to! 



- 1 



1 



:exp(-k'R), (32) 



where the normalization A/2 is defined in Eq. (|A4|l . As 
(i approaches there is a quantum phase transition be- 
tween the non-Abelian phase and Abelian phase. This 
transition is accompanied by closing of the gap and the 
Majorana bound state is no longer localized since k' — » 0. 

We briefly comment on the degeneracy splitting be- 
tween vortex zero modes in the ferromagnetic insula- 
tor/semiconductor /superconductor hybrid structure pro- 
posed by Sau et. al. [3l| which can be modeled by spin- 
1/2 fcrmions with Rashba spin-orbit coupling and s-wave 
pairing induced by the superconducting proximity effect. 
Since time-reversal symmetry is broken by the proximity- 
induced exchange splitting, this system belongs to the 
same symmetry class as spinless p x + ip y superconductor 
- class D. The connection between this hybrid structure 
and spinless p x + ip y can be made more explicit by the 
following argument: the single particle Hamiltonian after 
diagonalization yields two bands. Assuming a large band 
gap (which is actually determined by exchange field) , one 
can project the full Hamiltonian onto the lower band and 
then the effective Hamiltonian takes exactly the form of 
spinless p x + ip y superconductor , see, for example, the 



discussion in Ref. [32j. Although analytical expression 



for Majorana bound state in vortex core is not available, 
the solution behaves qualitatively similar to the one in 
spinless p x + ip y superconductor. Therefore, we expect 
that splitting should also resemble that of spinless p x +ip y 
superconductor . 



B. Splitting in TI/SC heterostructure 

In this section we discuss the case of vortex- vortex pair 
in TI/SC heterostructure. We assume both vortices have 
vorticity 1. Similar the case of p x + ip y superconductor, 
one can transform the surface integral over half plane S 
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to a line integral along its boundary 9S. Exploiting the 
explicit expressions for the zero mode solution, we arrive 
at 



/ d 2 r *l«BdG*+ - / d 2 r^H B dG*i 



2\[2v / dyxt( s )x;( s ) c °S¥'2, 



(33) 



is no splitting between two vortices with the same vortic- 
ity from this line of reasoning. As discussed below this 
fact actually holds beyond perturbation theory and the 
robustness of zero modes in the presence of chiral sym- 
metry is ensured by an index theorem. 



IV. ATIYAH-SINGER-TYPE INDEX THEOREM 



where s — \J {R/2) 2 + y 2 , cos ^2 = R/2s. 

First we consider the case with finite /i. There are two 
length scales: Fermi wavelength kp 1 = — and coherence 
length £ = . We evaluate the integral (f33|) in the limit 
where R is large compared to both kp 1 = v//i and £: 



£4 



4Afiv 



cos(kpR ■ 



0rfcF(l + 4£ 2 ) 1/4 ^Wl 



exp 



(34) 

where 2a = arctan(fc^£) and the normalization A3 is 
given by Eq. (|A7j) . One can notice that the splitting, in- 
cluding its sign, oscillates with the intervortex separation 
R when R is large. In the limit of large fi, say kp^ ^> 1, 
Eq. (|3~4"|) can be simplified to 



2A cos(k F R+ f) 



■exp - — 



(35) 



We now turn to the limit where fi is very close to Dirac 
point, i.e. /i — > 0, kp£, <C 1. We evaluate the integral for 
\ < R < kp 1 . 



E 4 



2[i 



3/2 



-r exp - — 



(36) 



where we have made use of asymptote of A/3 in the limit 
kp£ <C 1. Eq. (|31)j) implies that for fixed R the splitting 
vanishes as (i approaches Dirac point. Actually this fact 
can be easily seen from (1331) without calculating the in- 
tegral. Because at fi — either xt or xi vanishes, the 
splitting which is proportional to the product of Xi an d 
Xt is zero. The same result for splitting at fi = has also 
been obtained in Ref. [49j . 

We now show that vanishing of the splitting at fi = 
is a direct consequence of chiral symmetry. At fi = 
zero modes carry chirality which labels the eigenvalues 
of 7 5 . More specifically, wave function is an eigenstate 
of 7 5 : j 5 ^i — A^i- Consider an arbitrary perturbation 
represented by O to the ground state manifold expanded 
by these local zero modes. To leading order in pertur- 
bation theory its effect is determined by matrix element 
Oij = (^i\0\^j). Now assume that *i and *j have the 
same chirality (which means that vortices i and j have 
identical vorticity). If the perturbation O preserves chi- 
ral symmetry, i.e. {7 s , O} = 0, then 



(* 1 |{ 7 5 ,0}|* J > =2\(%\0\* j 



0. 



(37) 



Therefore matrix element | C| \I/ ^ ) vanishes identically. 
Tunneling obviously preserves chiral symmetry so there 



Index theorem provides an intelligent way of under- 
standing the topological stability of zero modes. It is 
well-known that one can relate the analytical index of an 
elliptic differential operator (Dirac operator) to the topo- 
logical index (winding number) of the background scalar 
field in 2D j46|) through the index theorem. Since BdG 
Hamiltonian for TI/SC system at /1 = can be presented 
as a Dirac operator (see Eq. p^)) ). we give a brief review of 
this index theorem, see also recent exposition in Ref. [73| . 
Specifically, the Hamiltonian for TI/SC heterostructure 
can be written as 



U D = n • V + T • n, 



(38) 



where n = (ReA, — ImA) field describes the non-trivial 
configuration of the superconducting order parameter. 
We assume the following boundary condition for n field: 



|n(r)| 



const as r 



(39) 



As mentioned above, this model Hamiltonian has 
particle-hole symmetry, time-reversal symmetry and chi- 
ral symmetry which is given by 7 s . It anticommutes with 
the Dirac Hamiltonian {j 5 ,Hd} = 0. Therefore, all zero 
modes ^0 of Hd are eigenstates of 7 5 . Since (7 5 ) 2 = 1 
eigenvalues of 7 s are ±1. We define ± chirality of zero 
modes as 7 5x &J = i'Sj- The analytical index of T-Lb is 
defined as 



ind"H/ 



(40) 



where n± are number of zero modes with ± chirality. 

The index theorem for the Hamiltonian Ho states that 
the analytical index is identical to the winding number 
of the background scalar field in the two-dimensional 
space : 



indHn = - — 
Ztt 



d 1 xe ab h a d i n b , 



(41) 



where n = n/|n|. According to the index theorem, the 
number of zero modes is determined by the topology of 
order parameter at infinity. The right hand side is en- 
sured to be an integer by the fact that the homotopy 
group tti(S 1 ) = Z. If we have a vortex in the system 
with vorticity I, the right hand side of (|4"Tj) evaluates ex- 
actly to /. Thus the index theorem implies that the Dirac 
Hamiltonian has at least / zero modes which a gree s with 
explicit solution obtained by Jackiw and Rossi 45] . This 
conclusion can be generalized to the case where multiple 
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vortices are present. In that case the right hand side is 
basically the sum of vorticities of all vortices. 

The index theorem (pTTj) requires chiral symmetry 
which is broken by presence of a finite chemical poten- 
tial / 0. Now we argue that when chiral symmetry 
is broken the Majorana zero modes admit a Z2 classi- 
fication corresponding to even-odd number of zero en- 
ergy solutions. Generally speaking, a small chiral sym- 
metry breaking term cause coupling between zero modes 
and split them away from zero energy. However, due to 
particle-hole symmetry, the number of zero modes that 
are split by chiral symmetry breaking term must be even. 
So the parity of the topological index is preserved in the 
generic case. This is consistent with an explicit solu- 
tions of zero mode in TI/SC heterostruture with finite 
chemical potential. Thus, we conclude that without chi- 
ral symmetry the Majorana zero modes bound to vortices 
are classified by Z2 corresponding to even or odd number 
of zero modes. 

Now we can fit our splitting calculation into the gen- 
eral picture set by index theorem. As being argued above, 
Majorana zero modes in spinless p x + ip y superconduc- 
tor is classified by Z2. When there are two vortices in 
the bulk, the topological index of order parameter is 2 
thus there is no zero mode and we find the splitting as 
expected. The same applies to two vortices in TI/SC 
heterostructure with (i / 0. However, as we have seen 
in the calculation the splitting vanishes for /1 = 0. This 
should not be surprising since according to index theo- 
rem, there should be at least two zero modes associated 
with total vorticity 2 which is the case for two vortices. 



A. Comparison with the splitting calculations in 
other systems. 

Recently numerical calculations of the degeneracy 
splitting have been performed for other systems support- 
ing non-Abelian Ising anyons [HTl I52I . l54j . In all these 
calculations it was found that the splitting has qualita- 
tively similar behavior - there is an exponential decay 
with the oscillating prefactor which stems from the spa- 
tial oscillations of Majorana bound states. In the case 
of Moore- Read quantum Hall state [52J , the splitting be- 
tween two quasiholes exponentially decays and oscillates 
with the magnetic length l c = since there only one 
length scale in the problem. The oscillatory behavior 
is also predicted for pair of vortex excitations in the B- 
phase of Kitaev's honeycomb lattice model in an external 
magnetic field |54j . 



V. COLLECTIVE STATES OF MANY-ANYON 
SYSTEM 

The microscopic calculations of the degeneracy split- 
ting for a pair of vortices are important for understand- 
ing the collective states of anyons arising on top of the 



non-Abelian parent state when many Majorana fermions 
(Ising anyons) are present [H, [73 - l77| . Essentially, the 
sign of the splitting favors certain fusion channel ( 1 or 
ip in the terminology of Ref . [lj ) when two vortices carry- 
ing Majorana fermions are brought together. These fu- 
sion channels correspond to having a fermion (^-channel 
when E+ < 0) or no fermion (1-channel when E + > 0) 
left upon fusing of two anyons. 

For pedagogical reason we start with the dilute anyon 
density limit assuming that the average distance between 
Majorana fermions is large compared with the coher- 
ence length £. In this regime, the many anyon state of 
the system will resemble gas of weakly bound pairs of 
anyons formed out of two anyons separated by the small- 
est distance. Because of the exponential dependence of 
the energy splitting the residual "interactions" with other 
anyons are exponentially smaller and can be ignored. In 
this scenario the parent state remains unchanged. 

When the density of anyons is increased so that the av- 
erage distance between them becomes of the order of the 
Majorana bound state decay length (coherence length £ 
in p-wave superconductors or magnetic length l c in Quan- 
tum Hall states) the system can form a non-trivial col- 
lective liquid (Wigner crystal of anyons or some other 
incompressible liquid state). This question has been in- 
vestigated in Refs. [6^, l78H80j . Although our approach 
used to calculate energy splitting breaks down in this 
regime and one should resort to numerical calculations 
for the magnitude of the splitting, we believe that qual- 
itative form of the splitting will remain the same. It 
is interesting to discuss the collective state that forms 
in this regime. Remarkably, it was shown in Ref. [69[ 
that depending on the fusion channel (i.e. sign of the 
splitting) the collective state of anyons may be Abelian 
or non-Abelian. This result was obtained assuming that 
the magnitude of the splitting is constant and the sign of 
the splitting is the same for all anyons (positive or nega- 
tive). However, because of the prefactor changing rapidly 
with the Fermi wave length we expect the magnitude of 
the splitting energy to be random realizing random bond 
Ising model discussed in Ref. [6^, [8l| . 

Finally, we mention that our calculations above and 
all existing studies of interacting many-anyon systems 
treat host vortices as classical objects with no internal dy- 
namics. This is a well-defined mathematical framework, 
which corresponds to the BCS mean-field approximation. 
In real superconductors, however, there are certainly cor- 
rections to it. The order parameter field, A(r, t), which 
describes a certain vortex configuration has a non-trivial 
dynamics and fluctuates in both space and time. At low 
temperatures, when the system is fully gapped, these 
fluctuation effects are suppressed in the bulk, but they 
are always significant in the vicinity of the vortex core, 
where the order parameter vanishes. This dynamics gives 
rise to an effective motion of a vortex as well as to the 
dynamics of its shape and the radial profile. The rel- 
evant length-scales of these effects certainly exceed the 
Fermi wave-length, which is the smallest length-scale in 
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the problem in most realistic systems. Even if the vor- 
tex is pinned, e.g. by disorder, its motion can be con- 
strained only up to a mean-free path or another rele- 
vant length-scale, which is still much larger than the 
Fermi wave-length for local superconductivity to exist. 
These considerations suggest that the intervortex separa- 
tion between quantum vortices has an intrinsic quantum 
uncertainty, which is expected to much exceed the in- 
verse Fermi wave- vector. This makes the question of the 
sign of Majorana mode coupling somewhat ill-defined in 
the fully quantum problem. Indeed we found the energy 
splitting to behave as SE(r) = \5Ea(r)\ cos {kpr + a), 
where |<5i?o(f)| is an exponentially small magnitude of 
coupling insensitive to any dynamics of r(t). The cosine- 
factor, which determines the sign, is however expected to 
be very much sensitive to quantum dynamics. To derive 
the actual microscopic model even in the simplest case of 
two non-Abelian anyons living in the cores of quantum 
vortices is a tremendously complicated problem, which 
requires a self-consistent treatment of the vortex order- 
parameter field and fermionic excitations beyond mean- 
field. However, one can argue that the outcome of such a 
treatment would be an effective theory where the e lkF " r ^ > 
factor that appears in Majorana interactions, should be 
replaced with a random quantum- fluctuating phase (c./., 
Ref. [82[), e l6< ^\ whose dynamics is governed by an ef- 
fective action of type, S[6] as J dr (9 — 6 ) 2 + c(d T 9) 2 . 
This generally resembles a gauge theory, but of an un- 
usual type, and at this stage it is unclear what collective 
many-anyon state such a theory may give rise to. 



VI. CONCLUSIONS 

In this paper, we address the problem of topological 
degeneracy lifting in topological superconductors char- 
acterized by the presence of Majorana zero-energy states 
bound to the vortex cores. We calculate analytically 
energy splitting of zero-energy modes due to the inter- 
vortex tunneling. We consider here canonical model of 
topological superconductor, spinless p x + ip y supercon- 
ductor, as well as the model of Dirac fermions coupled to 
superconducting scalar field. The latter is realized at the 
topological insulator/s-wave superconductor interface. 
In the case of spinless p x + ip y superconductor, we find 
that, in addition to the expected exponential decay, the 
splitting energy for a pair of vortices oscillates with 
distance in weak-coupling superconductor and these 
oscillations become over-damped as the magnitude of the 
chemical potential is decreased. In the second model, the 
splitting energy oscillates for finite chemical potential 
and vanishes at fi = 0. The vanishing of splitting 
energy is a consequence of an additional symmetry, the 
chiral symmetry, emerging in the model when chemical 
potential is exactly equal to zero. We show that this 
fact is not accidental but stems from the index theorem 
which relates the number of zero modes of the Dirac 



operator to the topological index of the order parameter. 
Finally, we discuss the implications of our results for 
many-anyon systems. 
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Appendix A: Normalization of Majorana bound 
state wave function 



In this appendix we present expressions for the normal- 
ization constants of Majorana bound state radial wave 
functions in Eqs.©,© and (p?Tj) . These constants are 
expressed in terms of hypergeometric functions. 

Normalization constant A/"i appearing in Eq.flT]) is de- 
fined as 



/>oo 

4nAf 2 / r&r Jl{kir)e 
Jo 



-2r/£ 



1, 



(Al) 



where fci = — A'^/v'p. Evaluation of the integral 

yields: 



A/? 



37r/cK 4 2 Ji(§,§; 3; 
with its asymptotes given by 



A/? 



< 1 



j| fci£»l 



(A2) 



(A3) 



Now we turn to A/2. Similarly, it's determined by 

POO 

AirNi / rdrI 2 (k 2 r)e- 2r/ Z = 1, (A4) 
Jo 



where fc 2 = \J A^jvp — 2m/i. Since ^ > 0, fc 2 £ is always 
smaller than 1. We find A/2 is given by 



A/? = 



3nk 2 e 2 ^i(|,|;3;fc|e 2 ) 



(A5) 



Finally, the normalization constant of wave function in 
(f!HT) can be calculated from 



-2r/S _ 



1, (A6) 



11 



which yields 



< 2 [82*1(5. 1; 1; -A 2 ) + 3A2 2 ^(§, |; 3; - A 2 )] 

(A7) 

with A = fj^/v. It has the following asymptotes: 



2? 
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